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ABSTRACT 

Context. The K2.5 dwarf HD 40307 has been reported to host three super-Earths. The system lacks massive planets and is therefore 
a potential candidate for having additional low-mass planetary companions. 

Aims. We re-derive Doppler measurements from public HARPS spectra of HD 40307 to confirm the significance of the reported 
signals using independent data analysis methods. We also investigate these measurements for additional low-amplitude signals. 
Methods. We used Bayesian analysis of our radial velocities to estimate the probability densities of different model parameters. We 
also estimated the relative probabilities of models with differing numbers of Keplerian signals and verified their significance using 
periodogram analyses. We investigated the relation of the detected signals with the chromospheric emission of the star. As previously 
reported for other objects, we found that radial velocity signals correlated with the S-index are strongly wavelength dependent. 
Results. We identify two additional clear signals with periods of 34 and 51 days, both corresponding to planet candidates with 
minimum masses a few times that of the Earth. An additional sixth candidate is initially found at a period of 320 days. However, this 
signal correlates strongly with the chromospheric emission from the star and is also strongly wavelength dependent. When analysing 
the red half of the spectra only, the five putative planetary signals are recovered together with a very significant periodicity at about 200 
days. This signal has a similar amplitude as the other new signals reported in the current work and corresponds to a planet candidate 
with M sin i ~ 7 M e (HD 40307 g). 

Conclusions. We show that Doppler measurements can be filtered for activity-induced signals if enough photons and a sufficient 
wavelength interval are available. If the signal corresponding to HD 40307 g is a genuine Doppler signal of planetary origin, this 
candidate planet might be capable of supporting liquid water on its surface according to the current definition of the liquid water 
habitable zone around a star and is not likely to suffer from tidal locking. Also, at an angular separation of ~ 46 mas, HD 40307 g 
would be a primary target for a future space-based direct-imaging mission. 

Key words. Methods: Statistical, Numerical - Techniques: Radial velocities - Stars: Individual: HD 40307 



1. Introduction 

Current high-precision spectrographs, such as the High 
Accuracy R adial Velocity Planet Searcher (HARPS; 
IMavor et all I2003I) a nd the High Resolution Echelle 
Spectrograph (HIRES; IVogt et al.L |1994|) . enable detections 
of low-mass planets orbiting nearby stars. During recent years, 
radial velocity (RV) planet searches have revealed several 
systems of super-E arths and/or Neptun e- mass planets around 
nearby stars (e.g. Mavor et all l2009allbl : ILovis et a! 1 120111 
iPeoe et all [201 lb iTuomiL 120121 

The system of three super-Earths orbiting HD 40307 has re- 
ceived much attention because the planets appear i n dynamically 
packe d orbits close to mean motion resonances (IMavor et al.L 
2009a). This has been used as an argument to suggest that low- 
mass planets may be found in highly compact multiple systems 
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that are still stable in long-term, e.g. a possibility of having ten 
planets with masses of 17 M m within a distance of 0.26 AU on 
stable orbits dFunk et all 1201 Oh . However, the physical nature of 
these comp anions as scaled-up v ersions of the Earths is not en- 
tirely clear (IB antes et al 1 12009al) . Their masses, between those 
of Earth and Neptune, suggest they are Neptune-like proto-gas 
giants that could not accumulate enough gas before it was blown 
away by the newly born star. On the other hand, recent tran- 
sit observations of hot super-Earths around bright nearby stars 
(ILeger et all 120091 iBatalha eTaill201 UlWinn et : alll2CTTlh indi- 
cate that a good fraction of these hot super-Earth mass objects 
can have rocky compositions. 

In this article we re-analyse the 345 HARPS spectra publicly 
available through the ESO archive using a newly developed soft- 
ware tool called HARPS -TE RRA (template-enhanced radial ve - 
locity re-analysis application; Anglada-Escude & Butletil2012h . 
Instead of the classic cross-correlation function method (CCF) 
implemented by the standard HARPS-ESO data reduction soft- 
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ware (HARPS-DRS), we derive Doppler measurements by least- 
squares matching of each observed spectrum to a high signal-to- 
noise ratio template built from the same observations. A descrip- 
tion of the method and th e imple mentation details are given in 
lAngl ada-Esc ude & Butle] d2012h . In addition to an increase in 
precision (especially for K and M dwarfs), this method allows us 
to perform additional analyses and tests beyond those enabled by 
the CCF data products provided by the HARPS-DRS. As an ex- 
ample, it allows us to re-obtain the RV measurements using only 
a restricted wavelength range. As we show in Section|5j this ca- 
pability can be instrumental in ruling out the planetary nature of 
prominent signals correlated with stellar activity. 

We rely on the Bayesian framework when estimating the or- 
bital parameters supported by the data, determining the signifi- 
cances of the signals, and the modelling of the noise in the mea- 
surements. In previous studies, radial velocities received within 
an interval of an hour or so have been commonly binned together 
in an attempt to reduce the noise caused by stellar- surface- 
related effects, i.e. stellar oscil lations and granu l ation, and other 
factors within this timescale (Dumusque etaTL l20Toh . In prin- 
ciple, this would enable the detections of planets smaller than 
roughly 5 M© with HARPS over a vari ety of orbital distances, 
even a t or near the stellar habitable zone (Dumusque e t al.Ll2010l 
[20TTah . In our approach, and instead of binning, we apply a self- 
consistent scheme to account for and quantify correlated noise 
in the Bayesian framework and use Bayesian model probabili- 
ties to show that a solution containing up to six planets is clearly 
favoured by the data, especially when the redmost part of the 
stellar spectrum is used in the RV analysis. Only the confluence 
of refinements in these data analysis methods (re-analysis of the 
spectra and Bayesian inference) allows the detection and verifi- 
cation of these low-amplitude signals. 

We start with a brief description of the stellar properties of 
HD 40307 (Section 2) and describe the statistical modelling 
of the observations and the data analysis techniques we used 
(Section 3). In Section 4 we describe the properties of the RV 
measurements and perform a detailed Bayesian analysis that 
identifies up to three new candidate signals. We discuss the stel- 
lar activity indicators and their possible correlations with the RV 
signals in Section 5. In this same section, we find that one of the 
candidates is spuriously induced by stellar activity by showing 
that the corresponding periodic signal (P ~ 320 days) is strongly 
wavelength dependent. When the RVs obtained on the redmost 
part of the spectrum are analysed (Section 6), the 320-day signal 
is replaced by a signal of a super-Earth-mass candidate with a 
200-day period with a minimum mass of about ~ 7M© orbiting 
within the liquid water habitable zone of HD 40307. The anal- 
ysis of the dynamical stability of the system (Section 7) shows 
that stable solutions compatible with the data are feasible and the 
potential habitability of the candidate at 200 days (HD 40307 g) 
is discussed in Section 8. We give some concluding remarks and 
discuss the prospects of future work in Section 9. 

2. Stellar properties of HD 40307 

We list the basic stellar properties of HD 40307 in TableQ] This 
K2.5 V star is a nearby dwarf with a Hipparcos parallax of 77.95 
+ 0.53 mas, which implies a distance of 12 .83 + 0.09 pc. It is 
somewhat smaller (M* = 0.77 + 0.05 M e ; ISousa et all l2008j) 
and le ss luminous (log L stai /L = -0.639 + 0.060: iGhezzi et all 
201 Oh than the Sun . The star is quiescent (log^ K < -4.99; 
Mavor et all l2009al) and rela tively metal-poor with [Fe/H] = - 
0.31+0.03 (ISousa et all [2008). It also lacks massive planetary 
companions, which makes it an ideal target for high-precision 



Table 1. Stellar properties of HD 40307. 



Parameter 


Estimate 


Reference 


Spectral Type 


K2.5 V 


Grav et al. (2006) 


l°g*HK 


-4.99 


Mavor et al. (2009a - ) 


7i [mas] 


HH 0£_i_n CQ 


van Leeuwen (2007) 


log Lsku/Lq 


-0.639+0.060 


Ghezzi et al. (2010) 




4.47+0.16 


Sousa et al. (2008) 


M star [M G ] 


0.77+0.05 


Sousa et al. (2008) 


Teff [K] 


4956+50 


Ghezzi et al. (2010) 


[Fe/H] 


-0.31 ±0.03 


Sousa et al. (2008) 


V sin i [kms~'] 


<1 


Mavor et al. (2009a - ) 


P mt [days] 


-48 


Mavor et al. (2009a - ) 


Age [Gyr] 


-4.5 


Barnes (2007) 



RV surveys aiming at finding lo w-mass planets. According to 
the calibration of lBarnesI (120071) . HD 40307 likely has an age 
similar to that of the Sun (~ 4.5 Gyr). 

3. Statistical analyses 

3.1. Statistical models 

We modelled the HARPS RVs using a statistical model with a 
moving average (MA) term and two additional Gaussian white 
noise components consisting of two independent random vari- 
ables. 

The choice of an MA approach instead of binning is based on 
the results of lTuomi et all d2012bl) and accounts for the fact that 
uncertainties of subsequent measurements likely correlate with 
one another at time-scales of an hour in an unknown manner. We 
limit our analysis to MA models of third order (MA(3) models) 
because higher order choices did not improve the noise model 
significantly. Effectively, the MA(3) component in our noise 
model corresponds to binning. However, unlike when binning 
measurements and artificially decreasing the size of the data set, 
this approach better preserves information on possible signals in 
the data. The two Gaussian components of the noise model are 
the estimated instrument noise with zero-mean and known vari- 
ance (nominal uncertainties in the RVs) and another with zero- 
mean but unknown variance corresponding to all excess noise in 
the data. The latter contains the white noise component of the 
stellar surface, usually referred to as stellar "jitter", and any ad- 
ditional instrumental systematic effects not accounted for in the 
nominal uncertainties. Kep lerian signals and w hite noise com- 
ponent were modelled as in Tuo mi et all d201 lh . 

In mathematical terms, an MA(p) model is implemented on 
measurement m-, as 

p 

mi = r k {ti) + y + e-, + ^ (f>j[m H - r k {t^j) - y\ exp (t H - ?,■), (1) 

j=i 

where rjt(f,-) is the superposition of k Keplerian signals at epoch 
f; and y is the reference velocity. The random variable e, is 
the Gaussian white noise component of the noise model with 
zero-mean and variance cr 2 - cr 2 + <t\, where cr? is the (fixed) 
nominal uncertainty of the ith measurements and cr 2 is a free 
parameter describing the magnitude of the jitter component. 
Finally, the free parameters of the MA(p) model are denoted as 
(f>j,j = 1 , p - they describe the amount of correlation between 
the noise of the z'th and i - jth measurements. The exponential 
term in Equation <[TJ ensures that correlations in the noise are 
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modelled on the correct time-scale. Specifically, using hours as 
units of time, the exponential term (that is always < 1 because 
ti > U-j) vanishes in few hours. 

To demonstrate the impact of binning on this data set, 
Section 14.11 shows the analyses of binned RVs with the com- 
mon assumption that the excess noise is purely Gaussian. This 
corresponds to using the nightly average as the individual RV 
measurements and setting <pj = for all j in Eq. ([TJ. 

3.2. Bayesian analyses and detection thresholds 

To estimate the model parameters and, especially, their un- 
certainties as reliably as possible, we drew random samples 
from the parameter posterior densities using posterior sam- 
pling algorithm s . We used the adaptive Metropolis algorithm 
of lHaario et al.l d2001l) because it can be used to receive ro- 
bust samples from the posterior density of the parameter vec- 
tor when applied to models with m ultiple Keplerian signals (e.g. 
iTuomi etalt l20TTl iTuornl l2012h . This algorithm is simply a 
modified version of the famous Metr opolis-Hastings Markov 
chain Monte Ca rlo (MCMC) algorithm ([Metropolis et all 1 19531: 
lHastingsl 1 1 970b . which adapts the proposal density to the infor- 
mation gathered from the posterior density. We performed sam- 
plings of models with k = 0, 1, ...7 Keplerian signals. 

The samples from the posterior densities were then used to 
perform comparisons of the different models . We used the one- 
block Metropolis-Hastin gs (OBMH) method dChib & Jeliazkovi 
l200ltlcivde et al., 2007) to calculate the relative posterior prob- 
abilities of models with differing numbers of Keplerian sig- 
nals (e.g.lTuomi & KotirantaL 120091: iTuomil 1201 it ITuomi et all 
201 1; Tuomi, 2012). We performed several samplings using dif- 
ferent initial values of the parameter vector and calculated the 
means and the corresponding deviations as measures of uncer- 
tainties of our Bayesian evidence numbers P(m\Mk), where m 
is the measurement vector and Ait denotes the model with k 
Keplerian signals. 

The prior probability densiti es in our ana lyses were essen- 
tially uniform densities. As in ITuomil (1201 2|) . we adopted the 
priors of the RV amplitude n(Kj) = 1/(0, a«y), reference veloc- 
ity n(y) - U(-aRv,aRv), and jitter n{ 07) = U(0,a RV ), where 
U(a, b) denotes a uniform density in the interval [a, b]. Since the 
observed peak-to-peak difference in the raw RVs is lower than 
10 m s _1 , the hyperparameteraRy was conservatively selected to 
have a value of 20 ms _1 . The priors of the longitude of peri- 
centre (a>) and the mean anomal y (Mp) were set t o U (0 , 2n), 
in accordance with the choice of iFord & Gregory! d2007l) and 
ITuomil d2012l) . We used the logarithm of the orbital period as 
a parameter of our model because, unlike the period as such, it 
is a scale-invariant parameter. The prior of this parameter was 
set uniform such that the two cut-off periodicities were T m - m and 
Tmax- These hyperparameters were selected as T„,i„ = 1 .0 days 
and T max = lOT^, because we did not expect to find signals 
with periods less than 1 day. Also, we did not limit the period 
space to the length of the baseline of the HARPS time series 
(Tpbs X because signals i n excess of that can be detected in RV 
data dTuomi et al. , 2009) and because there might be long-period 
signals apparent as a trend with or without curvature in the data 
set. 

Unlike in traditional Bayesian analyses of RV data, we did 
not use uniform prior densities for the orbital eccentricities. 
Instead, we used a semi-Gaussian as n(ei) oc Af(0, cr^) with the 
corresponding normalisation, where the hyperparameter cr e was 
chosen to have a value of 0.3. This value decreases the posterior 
probabilities of very high eccentricities in practice, but still en- 



ables them if they explain t he data better than lower ones (Tuomi, 
12012b ITuomi et alll2012al) . 

For the MA components <pj, we selected uniform priors as 
n((f)j) = U(—l, 1), for all j = 1,2,3. This choice was made to 
ensure that the MA model was stationary, i.e. time-shift invariant 
- a condition that is satisfied exactly when the values are in the 
interval [-1, 1]. 

Finally, we did not use equal prior probabilities for the mod- 
els with differing numbers of Keplerian signals. Instead, follow- 
ing |Tuom3 (Hn!), we set them as P(M k ) = 2P(Mt+i), which 
means that the model with k Keplerian signals was always twice 
as probable prior to the analyses than the model with k + 1 
Keplerian signals. While this choice makes our results more ro- 
bust in the sense that a posterior probability that exceeds our 
detection threshold is actually already underestimated with re- 
spect to equal prior probabilities, there is a physical motivation 
as well. We expect that the dynamical interactions of planets in 
any given system make the existence of an additional planet less 
probable because there are fewer dynamically stable orbits. This 
also justifies the qualitative form of our prior probabilities for 
the eccentricities. 

Our criterion for a positive detection of k Keplerian signals 
is as follows. First, we require that the posterior probability of 
a fe-Keplerian model is at least 150 times greater than that of 
the k - 1 -Keplerian model dKass & Raftervl Il995t ITuomil l201ll 
2012; lTuomi et alll201 UlFeroz et alll201 lb . Second, we require 
that the radial velocity amplitudes of all signals are statistically 
significantly different from zero. Third, we also require that the 
periods of all signals are well-cons trained from ab ove and below. 
These criteria were also applied in TuomJ d2012l) . 

We describe the parameter posterior densities using three 
numbers, namely, the maximum a posteriori (MAP) estimate 
and the limits of the corresponding 99% Bay esian credibility sets 
(BCSs) or intervals in one dimension (e.g. ITuomi & KotirantaL 
I2OO9I) . 

3.3. Periodogram analysis 

As is traditionally the case when searching for pe riodic signals 
in time series, we used least-squares periodograms (Lombl ll97rt 
Scargle, 1982) to probe the next most significant periods left in 
the data. In particular, we use d the least-squares periodograms 
described in lCumm ing (2004), which adjust for a sine wave and 
an offset at each test period and plot each test period against the 
F-ratio statistic (or power) of the fit. While strong powers likely 
indicate the existence of a periodic signal (though strong powers 
may be caused by sampling-related features in the data as well), 
the lack of them does not nece ssarily mean t hat there are no sig- 
nificant periodicities left (e.g. iTuomlL 120 1 2l) . This is especially 
so in multi-Keplerian fits due to strong correlations and aliases 
between clearly detecte d signals and yet-undetected l ower am- 
plitude companions (e.g. Anglada-Escude. et alll2010l) . The rea- 
son is that residuals must necessarily be calculated with respect 
to a model that is assumed to be correct, which is clearly not the 
case when adding additional degrees of freedom, i.e. additional 
planetary signals, to the model. Therefore, determining the reli- 
ability of a new detection based on goodness-of-fit comparisons 
is prone to biases, which effe ctively reduce s the sensitivity and 
reliability of these detections (Tuomil [2012b - 

While periodograms of the residuals are very useful, they do 
not properly quantify the significance of the possibly remaining 
periodicities and, therefore, we used them as a secondary rather 
than a primary tool to assess the significance of new signals. The 
analytic false-alarm probability (FAP) thresholds as derived by 
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Cumming (120041) are provided in the figures as a reference and 
for illustrative purposes only. We used the same periodogram 
tools to assess the presence of periodicities in the time series of 
a few activity indicators. 

4. Analysis of the RV data 

The 345 measurements taken on 135 separate nights were ob- 
tained over a baseline of ~ 1900 days. In contrast to the discus- 
sion in M ayor et"a"fl (l2009ah . we could not confirm a long-period 
trend using our new RVs and we did not detect evidence for 
a trend in th e new CCF RVs obtained using t he HARPS -DRS. 
As shown in Anglada-Es cude & Butlerl (1201 2l) (Fig. 3), changes 
in the continuum flux accross each echelle order (also called 
blaze function) induce RV shifts of several ms _I if not prop- 
erly accounted f or. Th i s effec t was reported to affect HARPS 
measurements in iPepd ([2010) and appears to have been fixed 
by HARPS-DRS v3.5 based on the release note£] (issued on 
29 October, 2010). With respect to HD 40307 in particular, we 
found that when RVs were derived without blaze function cor- 
rection, a strong positive drift (~ 2 ma^'yr -1 ) was left in the 
Doppler time series. We speculate that the trend reported earlier 
may be caused by this blaze function variability. A notable fea- 
ture of the HD 40307 data set is that the epochs have a long gap 
of 638 days between 3055 - 3693 [JD-2450000]. This feature can 
effectively decrease the phase-coverage of the data on longer pe- 
riods and complicates the interpretation of periodograms due to 
severe aliases. 



4.1. Analysis of binned data 

In a first quicklook analysis, we worked with the nightly aver- 
ages of the radial velocities as obtained from HARPS-TERRA 
using the standard setup for K dwarfs. In this setup, all HARPS 
echelle apertures are used and a cubic polynomial is fitted to 
correct for the blaze function of each aperture. After this correc- 
tion, the weighted means v e of all RV measurements within each 
night e are calculated. As a result, we obtain internal uncertain- 
ties of the order of 0.3-0.4 ms 1 for the HARPS-TERRA RVs. 
Because of stellar and/or instrumental systematic errors, we ob- 
served that these individual uncertainties are not representative 
of the real scatter within most nights with five or more measure- 
ments. Using three of those nights we estimate that at least 0.6 
ms~ l must be added in quadrature to each individual uncertainty 
estimate. After this, the uncertainty of a given epoch is obtained 
as crj 1 = Yjj (o"f) _1 , where the sum is calculated over all expo- 
sures obtained during a given night . Finally, based on their long- 
term monitoring of inactive stars, iPepe et aTl d201 ll) inferred a 
noise level of 0.7 ms _1 to account for instrumental and stellar 
noise. After some tests, we found that adding 0.5 ms _1 in quadra- 
ture to the uncertainties of the nightly averages ensured that none 
of the epochs had uncertainties below the 0.7 ms _I level. The 
typical uncertainties of a single night derived this way were of 
the order of 0.8 ms -1 . These corrections are basically only well- 
educated guesses based on the prior experience with RV data and 
reported stability of the instrument. Therefore, one must be espe- 
cially careful not to over-interpret the results derived from them 
(e.g., powers in periodograms and significance of the signals). In 
the fully Bayesian approach, we treat the excess noise as a free 
parameter of the model, therefore the Bayesian estimates of the 
noise properties should in principle also be more reliable. 
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Period [days] 

Fig. 1. Least-squares periodograms of the binned HD 40307 ra- 
dial velocities for the residuals of the models with three (top) to 
six (bottom) periodic signals. The analytic 10%, 1%, and 0.1% 
FAPs are shown as horizontal lines. 



First, we re-analysed the nighly binned RVs to see whether 
we cou ld independently reproduce the results of iMavor et alj 
(2009a) when HARPS-TERRA measurements and our Bayesian 
methods wer e used. Assuming an unknown Gaussian noise pa- 
rameter (e.g. iTuomi et al.L 1201 lh in addition to the estimated 
measurement uncertainties, the posterior samplings and the cor- 
responding model probabilities easily revealed the three strong 
signals corresponding to periods of 4.3, 9.6, and 20.4 days. The 
residual periodogram of the three-Keplerian model revealed ad- 
ditional strong periodicities exceeding the 1% FAP level (Fig. 
Q] top panel) and we tested more complicated models with up 
to six Keplerian signals. Especially, we tested whether the addi- 
tional power present in the three-Keplerian model residuals (Fig. 
[1] top panel) at periods of 28.6, 34.8, 51.3, and 308 days, peaking 
above the 10% FAP level, are statistically significant by starting 
our MCMC samplings at nearby seed periods. 

The global four-Keplerian solution was found to correspond 
to the three previously known super-Earth signals and an addi- 
tional signal with an MAP period of 320 days. This period was 
bounded from above and below and its amplitude was strictly 
positive - in accordance with our detection criteria. The corre- 
sponding posterior probability of the four-Keplerian model was 
1.9xl0 5 times greater than that of a three-Keplerian one, mak- 
ing the 320-day signal significant. In addition to this signal, we 
could identify a 51 -day periodicity (Fig. [1] second panel) that 
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satisfied the detection criteria as well. Including this fifth sig- 
nal in the model further increased the model probabilitity by a 
factor of 6.6xl0 6 . We could furthermore identify a sixth signal 
with our six-Keplerian model, correponding to a period of 34.4 
days. However, even though the samplings converged well and 
the solution looked well-constrained, the six-Keplerian model 
was only five times more probable than the five-Keplerian one 
and would not be detected using our criteria in Section [3~2l Both 
of the two new significant signals had MAP estimates of their 
radial velocity amplitudes slightly lower than 1.0 ms" 1 - the sig- 
nals at 51 and 320 days had amplitudes of 0.70 [0.31, 1.09] ms" 1 
and 0.75 [0.38, 1.12] ms -1 , respectively, where the uncertain- 
ties are denoted using the intervals corresponding to the 99% 
BCSs. We note that the periodogram of sampling does not have 
strong powers at the periods we detect (see Fig. 1 in Mayor et al.L 
I2009ah . 

Given the uncertain nature of the signal at 34.4 days and 
the potential loss of information when using the nightly aver- 
aged RVs (artificial reduction of the number of measurements), 
we performed a complete Bayesian reanalysis of the full dataset 
(345 RVs), now including the aforementioned moving average 
approach to model the velocities. 

4.2. Analysis using all RV measurements 

The analyses of the unbinned data immedia t ely sh owed the three 
previously announced signals (M ayor et al.ll2009ah with periods 
of 4.3, 9.6, and 20.4 days. Modelling the data with the superpo- 
sition of k Keplerian signals and an MA(3) noise model plus the 
two Gaussian white noise components, our posterior samplings 
and periodogram analyses identified these signals very rapidly, 
enabling us to draw statistically representative samples from the 
corresponding parameter densities. 

The residual periodogram of this model (three Keplerians 
and MA(3) components of the noise removed) revealed some 
significant powers exceeding the 0.1% and 1 % FAP level at 320 
and 50.8 days, respectively (Fig. [2] top panel). Samplings of the 
parameter space of a four-Keplerian model indicated that the 
global solution contained the 320-day periodicity as the fourth 
signal and yielded a posterior probability for the four-Keplerian 
model roughly 1 .5 x 1 6 times higher than for the three-Keplerian 
one. The nature of this signal and its relation to the stellar activ- 
ity (SectionfS]) is discussed in Section lB31 

We continued by calculating the periodogram of the residu- 
als of our four-Keplerian model (Fig. [2] second panel) and ob- 
served a periodogram power that almost reached the 1 % FAP 
level at a period of 50.8 days. Including this fifth signal further 
increased the posterior probability of our model by a factor of 
6.4 xlO 5 . 

The residuals of the five-Keplerian model contained a peri- 
odicity at 34.7 days (Fig. [2] third panel) exceeding the 1% FAP 
level. The corresponding six-Keplerian model with this candi- 
date received the highest posterior probability - roughly 5.0x10 s 
times higher than that of the five-Keplerian model. Since the 
parameters of this sixth candidate were also well-constrained, 
we conclude that including the 34.7-day signal in the statistical 
model is fully justified by the data. 

We also attempted to sample the parameter space of a seven- 
Keplerian model but failed to find a clear probability maxi- 
mum for a seventh signal (see also the residual periodogram 
of the six-Keplerian model in Fig. |2l bottom panel). Although 
the periodicity of the seventh signal did not converge to a well- 
constrained probability maximum, all periodicities in the six- 
Keplerian model at 4.3, 9.6, 20.4, 34.7, 50.8, and 320 days were 
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Fig. 2. Least-squares periodograms of all 345 RVs of HD 40307 
for the residuals of the models with three (top) to six (bottom) 
periodic signals together with the analytic 10%, 1%, and 0.1% 
FAPs. 

still well-constrained, i.e. their radial velocity amplitudes were 
statistically distinguishable from zero and their periods had clear 
probability maxima. Because we were unable to constrain the 
parameters of a seventh signal using the posterior samplings, we 
cannot be sure whether the corresponding Markov chains had 
converged to the posterior density and cannot reliably calculate 
an estimate for the posterior probability of the seven-Keplerian 
model. Therefore we stopped looking for additional signals. 

From this analysis, we can state confidently that there are 
six significant periodicities in the HARPS-TERRA radial veloc- 
ities of HD 40307 when the whole spectral range of HARPS is 
used. As we show in the next section, one of them has the same 
period as the chromospheric activity indicator (S-index) and re- 
quires more detailed investigation. The analysis of all 345 RVs 
indicates that for these data binning appears to be a retrograde 
step in extracting periodic signals from the RV data. We infer 
that binning serves to alter measurement uncertainties and damp 
the significance levels of the periodicities in the data. 

5. Stellar activity 

We examined the time series of two activity indicators derived 
from the cross correlation function properties as provided by the 
HARPS-DRS. They are the bisector span (BIS) and the full- 
width at half-maximum (FWHM) of the CCF These indices 



5 



Mikko Tuomi et al.: Habitable-zone super-Earth candidate in a six-planet system around the K2.5V star HD 40307 



monitor different features of the average stellar line. Briefly, BIS 
is a measure of the stellar line asymmetry and should correlate 
with the RVs if the ob served offsets are caused by spots or plages 
rotating with the star dOueloz et al.Ll200Th . The FWHM is a mea- 
sure of the mean spectral line width. Its variability (when not in- 
strumental) is usually associated with changes in the convective 
patterns on the surface of the star. A thi rd index, the s o -calle d 
S-index in the Mount Wilson system dBaliunas et all 1 19951) . 
is automatically measured by HARPS-TERRA on the blaze- 
corrected one-dimensional spectra provided by the HARPS- 
DRS. The S-index is a measure of the flux of the Call H and 
K lines (A H = 3933.664 A and A K = 3 968.470 A, respectiv ely) 
relative to a locally defined continuum dLovis et all 1201 lah and 
is an indirect measurement of the total chromospheric activity of 
the star. For simplicity, the analysis of the activity indicators was 
performed throughout for the 135 nightly averaged values using 
sequential least-squares fitting of periodic signals that are each 
described by a sine-wave model (period, amplitude and phase). 



5.1. Analysis of the FWHM and BIS 

The BIS was remarkably stable (RMS ~ 0.5 ms _1 ) and the pe- 
riodogram of its time series did not show any significant pow- 
ers. Visual inspection of the time series for the FWHM already 
shows a very significant trend of 5.3 ms _1 yr _1 . The 345 mea- 
surements of the FWHM are listed in Table IA.2I A sinusoidal 
fit to this trend suggested a period of 5000 days or more (see 
top panel in Fig. [3]). After removing the trend, two more signals 
strongly show up in the residuals. The first one was found at 
23 days and had an analytic FAP of 0.005%. After fitting a si- 
nusoid to this signal and calculating the residuals, an extremely 
significant peak appeared at 1 170 days with an analytic FAP of 
0.002%. After including this in a model with three sinusoids, no 
additional signals could be seen in the periodogram of the resid- 
uals with analytic FAP estimates lower than 10% (Fig.[3j bottom 
panel). We also show the FWHM values together with the fitted 
periodic curves in Fig. [4] While the signals in the FWHM were 
significant, we did not clearly detect their counterparts in the 
RVs. Given that BIS does not show any obvious signals either, 
we suspect that the periodicities in the FWHM might be caused 
by instrumental effects, e.g. tiny changes of the focus inside the 
spectrograph, rather than intrinsic variability of the stellar lines. 
The instrumental origin would reconcile the absence of corre- 
spondent drifts in the BIS and in the RVs. A similar indication 
of drifts and sensitivi ty of the FWHM t o instrumental issues has 
been reported in e.g. lLovis et alJ (1201 lah . While this adds some 
caveats on the long-term stability of the HARPS instrumental 
profile (and therefore its long-term precision), unless it is found 
to have similar periods, we see no reason to suspect that any of 
the signals in the RVs are spuriously induced by changes in the 
FWHM (intrinsic or instrumental). 
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Fig. 3. Periodogram series of the signals detected in the FWHM, 
from most significant to less significant (top to bottom). 
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Fig. 4. Time series of the FWHM activity index. The solid black 
line represents the best fit to a model containing three sinusoids 
(periods of 5000, 23 and 1 170 days). 



5.2. Analysis of the S-index 

The S-index also shows strong coherent variablity. The full set 
of 345 measurements of the S-index are provided in Table IA.21 
Again, this analysis was performed for the nightly binned mea- 
surements using least-squares periodograms and the sequential 
inclusion of sinusoidal signals. The last two S-index measure- 
ments were well above the average and could not be repro- 
duced by any smooth function. Even if they are representative 
of a physical process, such outlying points cannot be easily 
modelled by a series of a few sinusoids because these would 



add many ambiguities to the interpretation of the results. When 
these two points were removed, the periodograms looked much 
cleaner and we could unambiguously identify clear periodicities. 
Therefore, the analysis discussed here is based on the first 133 
nightly averages only. As illustrated in the periodograms of Fig. 
[5] this analysis showed strong evidence for three signals in the S- 
index. From most significant to less significant, we found these 
to consist of a long-period quadratic trend (with a period ex- 
ceeding the data baseline), a signal with P~ 320 days and a third 
signal at 43 days. Fig. [6] shows the best fit to a parabolic trend 
together with the two periodicities at 320 and 43 days. As has 
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Fig. 5. Periodogram series of the signals detected in the S-index, 
from most significant to less significant (top to bottom). 



been reported for several other stars (Dumusqu e et al.L |201 lab , 
we find it likely that the long-period trend is the signature of 
the stellar cycle that is not yet fully covered by the baseline of 
the observations. This trend might also be related to the trend ob- 
served in the FWHM, perhaps through an increase of the average 
magnetic field. 

Using more epochs. ILo vis et all d201 lab reached the conclu- 
sion that the stellar magnetic cycle of HD 40307 is about 5000 
days - a result that is compatible with ours. However, the nature 
of the 320-day signal is less clear. This period seems to be too 
long to be caused by rotation because it is more than a factor of 
three longer than the long est currently know n rotation periods 
which are about 100 days (llrwin et al.Ll201 ll) . Long periods like 
this would be difficult to explain by angular momentum evolu- 
tion (IBarnes et all 120121) . Indeed, the third signal in the S-index 
(43 days) appears to be a much better match t o the expected stel- 
lar rotation of HD 40307 dMavor et al.Ll2009bl) . The S-index also 
shows a few other tentative periodicities at 26, 63, and 167 days 
(FAP ~ 2-5%). 

HD 40307 is a very quiet star (log ~ -4.99; 

Mavor et al., 2009a) and it is in the effective temperature range 
where the correlation between magnetic activity and RV jitter 
shows the smallest corr elation (~ 4300-5000 K, see Fig. 2 in 
iDumu saue et all 1201 lbl) . Therefore, apart from the two most ob- 
vious signals in the S-indices (long-period trend and 320-day pe- 
riod), the other tentative periodicities are unlikely to produce de- 
tectable effects on the RVs. Indeed, no RV counterpart was iden- 
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Fig. 6. Time-series of the S-index. The solid black represents the 
best fit to a model containing three sinusoids (periods of 4000, 
320, and 43 days). Long-period trend and two signals at 320 and 
43 days are clearly detected in the time-series of the S-index. 



tified at 43 days (or any of the potential fourth signals at peri- 
ods of 26, 63, or 167 days), not even at low significance. Several 
studies have demonstrated a clear correlation of the S-index with 
observed RV offsets. For example, Du musque et alJ (1201 lb ) per- 
formed a comprehensive analysis of simultaneous HARPS RVs 
and chromospheric emission measurements of G and K dwarfs, 
showing that the correlation between the chromospheric emis- 
sion and RV offsets can be clearly traced and quantified as a 
function of spectral tvpe. lPepe et alJ (1201 ll) used this correlation 
to detrend the RV measurements of HD 85512 (K6 V) and re- 
ported the detection of am sin i ~ 2M ffl candidate with an orbital 
period of ~ 58 days. More recently, lAnglada-Escude & Butler] 
d2012h showed that long-period signals correlated with Doppler 
signals can be strongly wavelength dependent even within the 
limited HARPS wavelength range. Since Keplerian signals can- 
not be wavelength depedent, this provides a new set of diag- 
nostics to distinguish activity-induced signals from Keplerian 
signals. In s tead o f using the detrending technique employed by 
Pepe et alJ (1201 ll) . we investigate in the next section if any of the 
signals are wavelength dependent. 



5.3. Wavelength-dependent signals 



Anglada-Escude & Butlei] d2012l) showed that there is a cor- 
relation between the S-indices and RV measurements of HD 
69830 (G8 V) and that an apparent long-period trend in the 
RVs completely disappeared when one obtains the velocities 
using the redmost half of the HA RPS spectra only. Following 
Anglada- Escude & Butlerl (1201 2l) . we call this effect chromatic 
jitter. Chromatic jitter can cause coherent RV variability (e.g., 
periodic blueshifts following the magnetic cycle) but can also 
behave randomly, effectively adding uncertainty and correlated 
noise to the measurements. At present, we do not sufficiently 
understand the underlying physics to be able to predict the ex- 
istence and nature of the wavelength dependence of activity- 
related signals. Despite this lack of detailed knowledge, the di- 
minished sensitivity of redder wavelengths to activity has been 
shown to be an efficient method in verifying or falsifying the 
existence of propo sed candidates. For example, CRIRES/VLT 
RV measurements (Huelamoet al., 2008; Figueiraet al., 2010) 
in the H-band (~ 1.5^m) were used to rule out the existence of 
a gas giant candidate around TW Hya that w as initially detected 
in the optical RVs by ISetiawan et alJ (120081) . Another example 
was the brown dwarf candidate around LP 944-20. This candi- 
date was first detected at optical wavelengths (2 kms -1 semi- 
amplitude) and then ruled out by the same group using nIR 
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Fig. 7. RMS as a function of the bluest echelle aperture used. 
The corresponding Abc are illustrated at the top of the panel. A 
shallow minimum is reached at aperture 26 (Abc = 433 nm), in- 
dicating activity induced signals/jitter at the bluest wavelengths. 



NIRSPEC/Keck observations dMartfn et all I2006I) . The chro- 
matic nature of activity-induced signals has also abundant qual- 
itative support from theo ry and simulations (e.g. iReiners et all 
l2010HBarnes et al.Ll201 lh ." Below we explore the wavelength de- 
pendence of the 320-day signal (and all the others) to investigate 

its (their) nature in more detail. 

As suggested by Anglada-Escude & Butlerl (I2012I) . the first 
qualitative way of assessing the presence of chromatic jitter con- 
sists of obtaining radial velocity measurements for all epochs 
using the stellar spectrum from 680 nm (the longest wavelength 
available at HARPS) to some blue cut-off wavelength Abc, and 
plot the RMS of the resulting RVs as a function of Abc- For a 
perfectly stable star, the RMS should decrease monotonically 
as more e chelle apertures are used to obtain the velocity data. 
However. lAnglada-Escude & Butlerl (I2012I) showed that even for 
quiet G and K dwarfs without reported planets, a minimum in 
this RMS was typically reached when Abc was chosen to be be- 
tween 450-550 nm. For M dwarfs, this cut-off was found to be 
closer to 600 nm. In Fig. [7] we show the RMS of the radial ve- 
locities of HD 40307 as a function of Abc and note that there ap- 
pears to be a shallow minimum at Abc =453 nm. Although this 
minimum is not very deep (no Keplerian solution was subtracted 
from the RVs), its existence already suggests the presence of 
wavelength-dependent noise that acts stronger in the blue part of 
the spectrum. 

The next test consists of assessing if a signal observed us- 
ing the full spectrum is also present when examining a sub- 
section of it. For simplicity, we restrict this test to the nightly 
binned measurements and only use periodogram tools. In agree- 
ment with !4.2l when the full spectrum was used, the fourth dom- 
inant signal was found at a period of ~ 320 days. Any activity- 
induced signals only present at the bluest wavelengths should 
lose their significance as one shifts Abc towards the red wave- 
lengths. As illustrated in Fig. [8] this behaviour is clearly found 
for the 320-day signal. Interestingly, when Abc is set to 430 nm 
(20th HARPS aperture), the 320-day peak becomes as high as a 
new one emerging around 120 days. When the Abc is set to 503 
nm (40th HARPS aperture), the 120-day signal clearly domi- 
nates the periodogram. It is noteworthy that the powers of both 
the 34- and 51 -day peaks remain mostly constant, which shows 
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Fig. 8. Residual periodogram of the three-planet solution as 
a function of the blue cut-off wavelength Abc- The top peri- 
odogram is obtained using the full-spectrum RVs. The bottom 
periodogram is obtained using the aperture 40 (505.0 nm) as 
Abc- All signals discussed in the text are marked with narrow 
vertical lines. The horizontal line represents the peak value of 
peak closer to 320 days. When Abc - 433 nm (central peri- 
odogram), the emerging 125-day peak is already as high as the 
activity-induced signal. The 120-day peak is an alias of the most 
likely period of 200 days for planet candidate d. 



that their significance is not seriously affected by the choice of 
the cut-off wavelength. 

As we show in the next section, the 120-day signal is a yearly 
alias of the most favoured period at 200 days, which also appears 
conspicuously in the bottom periodogram of Fig. [8] To further 
illustrate the wavelength dependence of the 320-day signal, we 
show the fitted RV amplitudes of the four candidate periodicities 
(periods of 34, 51, 200, and 320 days) as a function of Abc (Fig. 
[9). These amplitudes are obtained by simultaneous least-squares 
fitting of seven circular orbits to the seven periods of interest 
(periods are also allowed to adjust). While the planet candidates 
at 34 and 5 1 days show roughly constant amplitudes, the 320- 
day signal almost completely vanishes when Abc is redder than 
450 nm (K< 0.4 ms _I ). It is also worth noting that the 200-day 
signal increases its amplitude as the 320-day signal disappears. 
The 200-day signal is seriously affected by yearly aliases which, 
effectively, correlates its amplitude with the signal at 320 days. 
The 320-day period is so suppressed at red wavelengths that for 
any Abc redder than 450 nm, we were forced to keep the period 
fixed at 320 days to estimate its amplitude and to prevent it from 
converging to a very different periodicity. 

At this point, it is fair to ask why the long-term variability in 
the S-index (and in the FWHM) does not show in the RV data 
as well. We can show that it actually does, very prominently, 
when calculating the RVs using the bluemost part of the spec- 
trum only. To illustrate this, we obtained the RV measurements 
using the wavelength range from 380 nm to 440 nm (20 bluest 
echelle apertures). After subtracting the signals of the three plan- 
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Table 2. Overview of the significant periodicities found in two 
activity indicators 
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Fig. 9. Amplitudes of the four signals under discussion as a func- 
tion of the blue cut-off Abc- While all proposed candidates keep 
a constant (or slightly increase) the derived semi-amplitude, the 
320-day signal is almost nonexistent (the period actually be- 
comes unconstrained) when wavelengths redder than 480 nm are 
used to derived a seven-planet orbital solution. 
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Fig. 10. Residuals of the three-planet fit when only the bluemost 
part of the spectra were used to extract the Doppler measure- 
ments. A quadratic trend similar to the one detected in the S- 
index is the most obvious signal left in these residuals. 



ets reported by May or et"a"D d2009ah . a parabolic shape is clearly 
left in the RV residuals. Fig.[l0]illustrates how well the parabolic 
trend matches the long-period trend also observed in the S-index 
(compare Figs. l6l and [T0T>. No blue RV counterparts of the other 
tentative periodicities in the S-index were identified in the data. 

As a conclusion to this section, we identified the chromatic 
nature of the ~ 320-day signal. We note that the internal er- 
rors of the individual spectra, given in the third column of Table 
IA.2I have a median value of 0.4-0.5 ms _1 . This number is sig- 
nificantly smaller than the measured intranight scatter and also 
smaller than the best RMS we could receive from any six- or 
seven-planet fit. This means that systematic noise (stellar and/or 
instrumental) is a dominant source of uncertainty and that ac- 
curacy is not lost by ignoring the Doppler information at the 
bluest wavelengths. A careful optimisation of Abc might pro- 
vide even more significant detections but, given the unknown 
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nature of this chromatic jitter, choosing Abc — 490.0«m (37th 
HARPS aperture) appears to provide a fair compromise. We sus- 
pect that the observed behaviour of HD 40307 RVs (also HD 
69830, HD 85512, and several o ther nominally stable stars, see 
lAnglada-Escude & Butleit 1201 2|) is a manifestation of the same 
processes observed in magnetically active stars but occuring at 
lower RV amplitude levels. The correlation between the S-index 
and the occurrence of the RV period at 320 days allows us to es- 
tablish a very likely connection between this activity -related sig- 
nal and a spurious RV signal and rule out the Keplerian nature of 
the corresponding candidate. We are in the process of analysing 
the chromatic dependence of stellar jitter using a more repre- 
senta tive sample (i.e. the sample analysed in iDumusque et all 
1201 Tab , the results of which will be presented in a future pub- 
lication. Given the substantial changes in the potential signals 
in the system (e.g., appearance of the candidate at a period of 
125 or 200 days), we perform a complete orbital re-analysis of 
the redmost half of the spectrum only (490-680 nm) in the next 
section. 



6. Analysis of the radial velocities from the 490-680 
nm spectra 

First, we performed a rapid analysis of the RVs obtained us- 
ing the redmost part of the spectra using least-squares peri- 
odograms. As was found to be the case when analysing the RVs 
from full spectra, thr ee additional periodicities compared to the 
dMavor et all l2009al) solutions could be extracted from the data 
if all orbital eccentricities were fixed to zero. Two of the signals 
had the same periods (34 and 5 1 days) and the third one showed 
at ~ 200 days (and on its yearly alias at 125 days). Although a 
model with circular orbits always looks more appealing aesthet- 
ically, we admit that this assumption cannot be justified using a 
pure least-squares approach without the use of prior information 
and/or physical constrains (e.g., orbital dynamics). As before, 
and in order to obtain a robust detection of all the candidates, 
we performed a full Bayesian analysis of the 345 red RV mea- 
surements using the methods and models described in Section 
EH 

The parameter spaces of models with k = 0, ...,6 were 
easy to sample using the adaptive Metropolis algorithm and 
the Markov chains converged rapidly to th e posterior densities . 
Again, the three candidates presented by May or et"aT1 ((2009a) 
with periods of 4.3, 9.6, and 20.4 days were extremely signifi- 
cant and easy to find using our MCMC samplings. 

After removing the MA(3) component of the noise and the 
signals in the three-Keplerian model, the power spectrum of the 
residuals showed very strong peaks at periods of 50.8, 82.1, and 
128.3 days (Fig.QT| top panel) exceeding 0.1% FAP levels. Even 
though the highest peak corresponds to the 125-day period, the 
MCMC samplings showed that the 50.8-day periodicity is the 
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Fig. 11. As in Fig.[2]but for the unbinned radial velocities of the 
490-680 nm spectrum. 



global solution by providing higher posterior probabilities to the 
four-Keplerian model. We note that when calculating the resid- 
ual periodogram of the three-Keplerian model, it is necessary to 
assume a fixed 3-Keplerian solution (and noise model) for the 
data. In contrast, the MCMC automatically accounts for corre- 
lations by also adjusting the noise parameters simultaneously to 
find the most probable areas of the parameter space (not only 
minimising the least-squares sum). After several MCMC sam- 
plings starting in the vicinity of the different possible periods 
suggested by the periodogram, this 50.8-day signal was always 
found to correspond to the global four-Keplerian solution and 
made the four-Keplerian model 3.1 x 10 12 times more probable 
than the three-Keplerian one. 

The residual periodograms of four- and five-Keplerian mod- 
els had moderate powers at 34.6, 81, 125, and 197 days (Fig. 
[TTl panels two and three). We performed similar samplings of 
the parameter spaces of five- and six-Keplerian models using the 
favoured periodogram powers as starting points of the Markov 
chains. As a result, we found that the period space near the 197- 
and 34.6-day powers provided the global maxima and improved 
the posterior model probabilities by factors of 3.6 x 10 11 and 
2.5 X 10 9 when moving from k — 4 to k — 5 and from k = 5 to 
k = 6, respectively, where k is the number of Keplerian signals 
in the model. Removing all six signals and the MA(3) compo- 
nents of the noise from the data left the periodogram of the resid- 
uals without any considerable powers (Fig. [TT] bottom panel). 
Similarly as for the full spectrum, the period of the seventh sig- 



Table 3. Relative posterior probabilities of models with k = 
0, ...,6 Keplerian signals (Aik) and the periods {P s ) of the sig- 
nals added to the model when increasing the number of signals 
in the model by one. P(d\Mk) is the Bayesian evidence and RMS 
denotes the common root mean square value of the residuals for 
comparison. 



k 


P(M k \d) 


log P(d\M k ) 


RMS [ms- 1 ] 


P s [days] 





1.6xl0- 12 ' 


-765.92 ± 0.06 


2.20 




1 


8.5X10" 106 


-715.20 ± 0.04 


1.88 


9.6 


2 


7.7xl0~ 74 


-640.92 ±0.13 


1.50 


20.4 


3 


3.5X10" 34 


-548.90 ± 0.27 


1.16 


4.3 


4 


l.lxlO- 21 


-519.43+0.14 


1.05 


50.8 


5 


4.1xl0- 10 


-492.12 ±0.05 


0.99 


198 


6 


~ 1 


-469.40 ±0.12 


0.92 


34.7 



nal in a seven-Keplerian model did not converge to a clear prob- 
ability maximum, thus failing to satisfy our requirement that the 
period of a signal has to be well-constrained for a positive detec- 
tion. 

To summarise, the posterior probabilities (and the corre- 
sponding Bayesian evidences) of models with k = 0,1,2,... 
Keplerian signals were found to be heavily in favour of k = 6 
(Table [3]). This was the case although the selected prior proba- 
bilities penalise the models more heavily as the number of sig- 
nals increases. We do not show the results for a seven-Keplerian 
model in Table [3] because the posterior samplings did not con- 
verge to a clear maximum (nor several maxima) in the seven- 
Keplerian model and therefore we cannot be sure whether the 
OBMH estimate yielded trustworthy values for the marginal in- 
tegrals needed to assess the model probabilities. We plot the 
phase-folded MAP estimated orbits of the six Keplerian signals 
in Fig. IT21 after removing the MA(3) components. 

One of our detection criteria for Keplerian signals in the ra- 
dial velocities is that their amplitudes should be statistically sig- 
nificantly different from zero. The MAP estimates and the cor- 
responding 99% BCSs of the six-Keplerian model parameters 
are shown in Table [4] The BCSs of the RV amplitudes clearly 
indicate that all six amplitude parameters indeed differed signif- 
icantly from zero in accordance with the model probabilities that 
imply the existence of six Keplerian signals in the data. This was 
clearly the case despite the 34.7- and 198-day signals having am- 
plitudes with MAP estimates below 1.0 ms -1 . Also, the periods 
of all signals were well-constrained and had clear MAP values 
and close-Gaussian densities (Table|4j Fig.fTSl. The distributions 
of orbital eccentricities are also shown in Fig.Q~3] In agreement 
with the quicklook detection using least-squares periodograms 
assuming circular orbits, these distributions demonstrate that all 
signals had eccentricities consistent with close-circular orbits. 

We also analysed the RVs derived from two independent 
blocks of the red part of the spectrum to assess if any of the 
signals was only present at certain wavelength ranges. The first 
block comprised the range between 465 and 553 nm (from now 
on Redl) and the second one covered the range between 553 to 
690 nm (Red2). Redl contains the richest part of a HARPS spec- 
trum of a mid-K dwarf in terms of Doppler information (higher 
signal-to-noise, but also higher density of lines). The mean un- 
certainty of Redl is approximately 0.6 ms _1 . As a result, we 
found that the same six signals were directly spotted even with 
periodograms only. The amount of Doppler information in Red2 
is lower and nominal uncertainties were found to be approxi- 
mately 1.0 ms _1 on average. Using sequential periodogram sub- 
traction of the signals we could recover the same six candidates 
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Fig. 12. Phase-folded MAP Keplerian signals of the six-planet solution with the other five signals and the MA(3) components 
removed from each panel. 



except that the 125-day alias was preferred against the 200-day 
signal for HD 40307 g. However, the Bayesian samplings in- 
cluding the complete model for the uncertainties still converged 
to the same six-planet solution (P g ~ 200 days). Given that there 
is no a priori reason to prefer one or the other (125 versus 200 
days, no activity signals appear to be related to any of them) 
and because the MCMC samplings insist on favouring the 200- 
day period, we conclude that the solution in Table [4] is the one 
preferred by the data. The significances and model probabilities 
using Redl and Red2 are obviously lower than when using the 
whole red half of the spectrum. Yet, with the exception of the 34- 
day candidate in Red2 (chains converged to the same solution, 
but model probability ratio was close to the 1:150 threshold), all 
signals could be confidently identified even using two subsets of 
the HARPS data (each contains 1/3 of the full HARPS wave- 
length coverage). These tests give us further confidence in the 
reality of the reported candidate signals. 

As a final consistency check, we calculated the model prob- 
abilities using different information criteria and compared them 
to the Bayesian model probabilities derived using our approxi- 
mation for the marginal integrals. First, we estimated the rela- 
tive posterior prob abilities using the Akaike information crite - 
rion (AIC; see e.g. iLiddlel 120071 iBurnham & AndersonL 120021) . 
which is known to provide only rough approximations for the 
marginal integrals. The AIC (actually its small-sample version, 
AICc) yielded posterior probabilities for the different models 
that did not differ by more than a factor of 5 from the proba- 
bilities received using the OBMH estimate (Table [3). Using the 
AIC in model selection yielded the same qualitative results the 
OBMH estimate did. A common criticism of the AIC is that does 
not penalise the mo re complex models as much as it should in all 
applications (e.g. [Burnham & An dersonL 120021 and references 
therein). To obtain a more reliable comparison, we repeated the 
model selection procedure using the Bayesian information crite- 
rion (BIC), som etimes called the Schwarz information criterion 
(ISchwarzlll978h . which is known to penalise more complicated 
models, e.g. models with more Keplerian signals in this particu- 



lar case, more heavily than the AIC. Even according to the more 
restictive BIC, the six-Keplerian model had the highest posterior 
probability and was 1.3xl0 5 times more probable than the five- 
Keplerian model. These results provide further confidence in the 
conclusion that six periodic signals explain the variations in the 
data. 



6. 1 . Noise properties 

In addition to ruling out wavelength-dependent signals, the 
Bayesian analysis of the RVs at different wavelength ranges 
allows the comparison of the noise properties. As described 
in Section I3.2I we modelled the radial velocity noise using a 
third-order MA model (three free parameters) together with a 
Gaussian white noise component (one free parameter). The pos- 
terior distributions and MAP values for these parameters quan- 
tify the significance and the magnitude of each component of the 
noise (white vs correlated). 

The MA parameters <f> j = 1, 2, 3, which describe the corre- 
lation between the error of the i-th and (i-j)th measurement, are 
clearly positive for the RVs from full spectra (Tabled. The noise 
of the measurements appears to correlate positively the strongest 
with the noise of the previous measurement and the effect de- 
creases for the measurements obtained before that. However, this 
correlation decreases for the RVs obtained from the redmost part 
of the spectra (490-680 nm). According to our results, this de- 
crease is balanced by a corresponding increase in the magnitude 
of the Gaussian white noise component (Tabled. Therefore, the 
noise in the radial velocities becomes "whiter" when excluding 
the bluest half of the spectra from the analyses. This result means 
that in addition to removing long-term activity-related signals, 
using redder wavelengths also helps reducing correlations in the 
noise that might cause biases to any estimates of the Keplerian 
signals if not accounted for. 

We note that the six signals in the velocities were rather in- 
dependent of the exact MA model used. We obtained consistent 
solutions with MA(2) and MA(4) models as well, though these 
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Table 4. Six-planet solution of HD 40307 radial velocities. MAP estimates of the parameters and their 99% BCSs. Parameter ed yn 
denotes the estimates for eccentricity after excluding solutions corresponding to dynamically instability in 10 6 year time scales. 



Parameter 


HD 40307 b 


HD 40307 c 


HD 40307 d 


P [days] 


4.3123 [4.3111,4.3134] 


9.6184 [9.6135, 9.6234] 


20.432 [20.408, 20.454] 


e 


0.20 [0.04, 0.34] 


0.06 [0, 0.17] 


0.07 [0, 0.18] 




0.16 [0.04, 0.25] 


0.05 [0, 0.13] 


0.05 [0, 0.13] 


K [ms^ 1 ] 


1.94 [1.67, 2.25] 


2.45 [2.17, 2.75] 


2.75 [2.40, 3.10] 


a> [rad] 


3.4 [2.4, 4.2] 


4.1 [-] 


0.3 [-] 


M [rad] 


4.3 [2.5, 5.9] 


5.3 [-] 


5.6 [-] 


m p sin i [M e ] 


4.0 [3.3, 4.8] 


6.6 [5.6, 7.7] 


9.5 [8.0, 11.2] 


a [AU] 


0.0468 [0.0445, 0.0492] 


0.0799 [0.0759, 0.0839] 


0.1321 [0.1255,0.1387] 




HD 40307 e 


HD 40307 f 


HD 40307 g 


P [days] 


34.62 [34.42, 34.83] 


51.76 [51.30, 52.26] 


197.8 [188.8, 203.5] 


e 


0.15 [0, 0.28] 


0.02 [0, 0.22] 


0.29 [0, 0.60] 


^dyn 


0.06 [0, 0.18] 


0.03 [0, 0.13] 


0.22 [0, 0.45] 


K [ms- 1 ] 


0.84 [0.53, 1.16] 


1.09 [0.77, 1.37] 


0.95 [0.65, 1.27] 


a> [rad] 


5.3 [-] 


6.2 [-] 


1.6 [-] 


M [rad] 


6.2 [-] 


0.5 [-] 


5.6 [-] 


trip sin i [M e ] 


3.5 [2.1,4.9] 


5.2 [3.6, 6.7] 


7.1 [4.5,9.7] 


a [AU] 


0.1886 [0.1782, 0.1969] 


0.247 [0.233, 0.258] 


0.600 [0.567, 0.634] 


y [ms~'] 


-0.47 [-0.71,-0.25] 






(Tj [ms _1 ] 


0.73 [0.59, 0.87] 







Table 5. Noise parameters obtained when using the RVs from 
the full spectra and the red part (490-680 nm) of the spectra only 
as denoted using the MAP estimates and the 99% BCSs. 



Parameter 


Full 


Red 


<Pi 

<p2 
<p3 

(Tj [ms~'] 


0.54 [0.28, 0.77] 
0.23 [-0.13, 0.58] 
0.28 [-0.06, 0.66] 
0.53 [0.41, 0.65] 


0.27 [0.03, 0.50] 
0.10 [-0.25,0.42] 
0.25 [-0.07, 0.54] 
0.73 [0.59, 0.87] 



noise models were found to have slightly lower posterior prob- 
abilities. This indicates, as can be seen in Table [5] that the MA 
components MA(p) with p > 1 do not affect the solution much 
because their 99% BCSs imply that they are consistent with zero. 

7. Orbital stability 

We investigated the long-term stability of the planetary sys- 
tem corresponding to our six-Keplerian solution (Table @). This 
analysis was performed in three steps. First, we checked the 
Lagrange stability criteria, which provide a first qualitative as- 
sessement of the allowed orbital configurations. Second, we nu- 
merically integrated orbits directly drawn from the MCMC sam- 
plings and assessed the fraction corresponding to stable orbits 
at the Myr time scale. As a last test, we analysed the stability 
of the six-planet configuration by testing different eccentricities 
and arguments of the node of each planet. 

7. 1 . Lagrange stability criteria 

As a rough initial test, we checked whether the subsequent plan- 
ets in the system satisfied the ap proximate Lagrange s t ability cri- 
terion expressed analytically in Barnes & Greenberg (2006) and 
applied for instance in iTuomil d2012h . We plotted the resulting 
areas of instability in Fig. [14]- if the i - 1th or i + 1th planet has 
orbital parameters within the shaded areas surrounding the /th 
planet, the system is likely unstable over the long term. While 
this criterion is but a rough approximation of reality because it 
does not take the stabilising and/or destabilising effects of or- 
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Fig. 14. Approximate Lagrange stability thresholds (shaded ar- 
eas) between each two planets (red circles) and the MAP orbital 
parameters of the six-planet solution (Table©. 



bital resonances into account and is only valid for two planets at 
a time, it provides simple guidelines on the prospects of stability 
in the system. 

According to the Lagrange stability criterion and the cor- 
responding analytical thresholds in Fig. [14] the six planet can- 
didates orbiting HD 40307 form a stable system if their or- 
bital eccentricities are close to or below their MAP estimates. 
Alternatively, it can be stated that if the planets exist and have 
masses approximately equal to the minimum masses in Table[4] 
they are all likely on close-circular orbits. Especially, the outer 
companion is unlikely to have orbtital eccentricity in excess of 
0.4 (Fig. m. 
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Fig. 13. Distributions estimating the posterior densities of orbital periods (P x ), eccentricities (e x ), and radial velocity amplitudes 
(K x ) of the six Keplerian signals. The dashed lines show the prior densities for comparison for the orbital eccentricities. The solid 
curves are Gaussian densities with the same mean (fi) and variance (<x 2 ) as the parameter distributions. Additional statistics, mode, 
skewness (yu 3 ), and kurtosis (/i 4 ) of the distributions are also shown. 



7.2. Numerical integration of MCMC samples 

To further test the dynamical stability of the system, we inte- 
grated the planetary orbits for 10 6 years, whic h is appropriate to 
exclu de 'almost all' unstable configurations (iBarnes & OuinnL 
l200l. In our integration s, we used the Bulirsch-Stoer integrator 
(Bulirsch & Stoer, 1966) that has been used before when inte- 



grating planetary orbits (e.g. iTuomi et al.L 120091 : iTuomil 1201 lb . 
We found that most of the unstable configurations corresponded 
to one (or more) of the eccentricities of the companions being 
higher than its MAP estimate and the system turned out to be 
unstable in 10 4 - 10 5 years. Specifically, we drew a sample from 
the posterior density of the orbital parameters and used each vec- 
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tor in this sample as an initial state of orbital integrations. We 
increased the sample size such that there were 100 initial states 
that could not be shown to correspond to unstable orbital config- 
urations. This subsample was found to correspond to 4% of the 
total sample. 

In Fig. [15] we plot the posterior densities of the planetary 
eccentricities together with the densities obtained from the sub- 
sample of solutions that were not found to be unstable. Fig. [15] 
states essentially the same thing as Fig. [14] That is, it shows 
that the stable configurations mostly correspond to orbital ec- 
centricities closer to zero than those derived from the Bayesian 
MCMC analysis. The corresponding updated estimates for or- 
bital eccentricities are shown in Tableland denoted using ed yD . 
While this experiment shows that a non-negligible fraction of 
the solutions allowed by the data corresponds to stable config- 
urations, we aimed to explore the source of the dynamical in- 
stabilities and the role of potential dynamical resonances in the 
system as well. 

7.3. Sources of (in)stability 

Owing to the dense packing of the planets in the HD 40307 sys- 
tem, one still has to rely on systematic numerical integrations 
to correctly investigate any aspects of its stability. For example, 
from a direct numerical integration of the MAP six-planet solu- 
tion given in Table [4] one quickly finds that the orbital evolution 
leads to collisional trajectories within a timescale of 10 3 years. 

As was already shown above, it is likely that the planets 
have close-circular orbits. This configuration is also supported 
by the confidence intervals given for the individual eccentrici- 
ties in Table [4] which essentially state that circular orbits cannot 
be ruled out by the data. Therefore, to slightly simplify the dy- 
namical analysis, we performed a new fit to the observations, for 
which we allowed only circular orbits. For this nominal solution 
(S 1 hereafter) we analysed the impact of changing the initial ec- 
centricity e and initial mean anomaly M of individual planets on 
the stability of the whole system. 

Starting from SI we changed only the initial e or M of 
an individual planet. We then integrated the obtained sys- 
tem for 1000 years usi ng the symplectic SABA2 integrator 
dLaskar & Robutelll2001l) with a time step t — 0.1 days. The sta- 
bility of t he system wa s afterwards estimated using a frequency 
analysis (lLaskarl[l993l) . For this we divided the solution into two 
parts and computed the mean motions n for all planets for each 
part separately. To obtain the measure of stability (D), we cal- 
culated the maximum relative change in the difference between 
both mean motions as 



D = max 



(2) 



For a stable system, the individual mean motions will not vary 
with time, i.e. n\ * n l 2 . The results are shown in Fig. [16] Here 
"red" denotes strongly chaotic orbits, while "black" corresponds 
to very regular ones. 

The left panel of Fig. [TrSlconfirms that the stability of the sys- 
tem is crucially dependent on the individual eccentricities of the 
planets. Only for initial eccentricities very close to zero regular 
motion can be found (D ~ 10~ 9 ). Using the right panel of Fig. [16] 
it is also possible to identify the reasons for dynamical stability 
in HD 40307. While the system is not sensitive at all for exam- 
ple to the initial mean anomaly of planet g, it crucially depends 
on M e and Mj. In each column of these planets one finds in Fig. 
[16] initial conditions that lead to very regular motion (D ~ 10~ 9 



for e.g. Mf(to) » 310°) and also conditions that lead to chaotic 
motion (D ~ 10~ 2 for e.g. M f (t ) » 267°). When analysing the 
results of the respective integrations, we found that the orbital 
periods of candidates e and f need to have a ratio of 3:2 to main- 
tain stable motion. This resonant configuration helps to avoid 
close encounters between the planets due to the libration of the 
resonance angle 9 = -2A e + 3/1/ - o> e . We note that this config- 
uration is also maintained over long time intervals. Integrations 
up to 1 Myr using these initial conditions did not show any signs 
of chaotic motion. 

In contrast, for highly unstable initial conditions the angle 
9 circulates, which leads very soon to collisions between some 
planets or the ejection of one of them. 

In conclusion to the discussion on the dynamics of the sys- 
tem we find that - despite the presence of many planet candidates 
in the system - there exists a substantial number of dynamically 
stable orbits allowed by the data. In particular, orbits with ec- 
centricities close to zero are dynamically favoured due to the 
existence of the stabilising 3:2 mean motion resonance between 
candidates e and f. We note that we have only tested solutions 
based on the minimum masses of the candidates. It is possible 
that increasing the masses of all planets by a small factor (e.g., x 
1 .2) makes the system unstable and could therefore impose inter- 
esting upper limits to the masses. Additional RV measurements 
are necessary to better constrain the eccentricities of the several 
candidates and derive more robust conclusions on the dynamical 
features of the system. 

8. Habitability and prospective follow-up 

The so-called liquid water habitable zone (HZ) is the area around 
a star where an Earth-like planet could support liquid water on its 
surface. For HD 40307, this HZ lies between 0.43 and 0.85 AU 
(s ee Fig.fTTb. This in terval was computed using the recipes given 
in Selsis et al. (2007) and is a function of stellar luminosity and 
effective temperatu re. We used the lum inosity and temperature 
estimates given by iGhezzi et~aT ] (120101) . With an orbital semi- 
major axis of 0.60 AU, the planet candidate HD 40307 g receives 
~ 62% of the radiation the Earth receives from the Sun and lies 
safely within the HZ of its stellar host. Even though the radiation 
is somewhat low compared to that received by the Earth, we note 
that the Earth lies actually reasonably close to the inner boundary 
of the Sun's HZ (Fig. E). 

As is the case with the previously reported planets in the HZs 
of nearby stars, additional observational information is neces- 
sary - on top of confirming the planetary nature of the three new 
signals we report in this work - to decide if HD 40307 g indeed 
can support water on its surface. In the meantime, and given that 
it would be a rather massive obje ct for a telluric planet, detailed 



climatic ( e.g.lBarnes et all 20121) and planet ary interior Simula- 



_ „ intl nli 

tions (e.g. lKorenagal 2010HStein eLaU |2011) can be used to esti- 
mate its suitability for supporting liquid water and, perhaps, life. 
Unlike the other previously reported candidate habitable planets 
around nearby M dwarfs (e.g., GJ 581 and GJ 667C), HD 40307 
g is located farther away from the centr al star and, like the h ab- 
itable zone planet candidate Kepler 22b (Boruc kTet al.Ll2012l). it 
should not suffer fro m tidal locking either (Kasting et al., 1993; 
lBarnesetal.ll2009bl) . 

Given its relatively long-period orbit, the transit proba- 
bility of planet candidat e HD 40307 g is very low (~0.6% 
ICharbonneau el alll2007h . Also, previous attempts to detect the 
transits of the inne r candidate HP 403 07 b (P « 4.3 days) have 
been unsuccessful (Gillon et aTl l2010l) . Therefore, if coplanarity 
is a common feature of planetary systems, the fact that HP 
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Fig. 15. Distributions estimating the orbital eccentricities given measurements alone (grey histogram; mean fi m and standard devi- 
ation <r m ), showing the values the eccentricities had during orbital integrations of 10 6 years that could not be shown to be unstable 
(white histogram; mean jij and standard deviation era), and their combination (red histogram). 
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Fig. 16. Stability analysis of the circular six-planet solution (SI). The panels show the influence of changing the initial eccentricity 
(left) and initial mean anomaly (right) on the stability of the complete planetary system. "Red" corresponds to highly unstable orbits 
while "black" denotes regular motion. Each of the 100 initial conditions per planet were integrated numerically for 1 kyr, for which 
we then computed the stability criterion D (Eq. [3J. 



40307 b is not transiting its star bodes ill for the transit chances 
of any outer companion. HD 40307 is a nearby star (13 pc) com- 
pared to the typical target of the Kepler mission (e.g., Kepler 22 
is located at ~ 200 pc). Because of this, the star-planet angu- 
lar separation can be as large as ~ 46 mas. Space-based dire ct- 
imaging missi ons such as Darwin /ESA (ICockell et all l2009t) or 
NASA/TPF (Lawso rTet all 12006) aim to reach the sensitivity to 
detect Earth-sized planets at an inner working angle of ~ 40 mas. 
Therefore, HD 40307 g would be the first habitable-zone candi- 
date in the super-Earth mass regime that could be targeted by the 
planned direct imaging observatories. 



9. Discussion 

After examining the wavelength dependence of the signals in the 
RVs of HD 40307 and their relation to the chromospheric ac- 
tivity indicator (S-index), we easi ly detect the t hree su per-Earth 
candidates previously reported by May or et"aD d2009al) using in- 



dependent data analysis methods. In addition to these, we also 
confidently detect three new periodic signals in the RVs of this 
K2.5 V dwarf star. Our analyses of the HARPS RVs indicates 
the existence of a system with six super-Earth candidates around 
this star (Table|3j Fig.fTTTi and our dynamical analyses show that 
long-term stable orbits are allowed within our credibility sets 
of the parameters. Our results make HD 40307 one of the few 
known pla netary systems with more than fiv e planet candidates: 
HD 10180 (lLoyis et al.Ll201 la[lTuomill2012h . the Kepler- 1 1 six- 
planet system (ILissauer et all 1201 ll) . and the solar system. 

A handful of potentiall y habitable planet s have been re- 
p orted to date. GJ 581 d dMavor et all l2009bl) and HD 88512 
b dPepe et alll201 ll) orbit at the extended HZ where thick clouds 
of CO2 or water would be require d to sustain wet atm ospheres. 
As is th e case with Kepler-22 b dBorucki et all 120121) and GJ 
667C c (lAnglada-Escude et all 120121) . HD 40307 g orbits well 
within the more classic and conservative definition of the HZ. 
With a radius of 2.4 /? ffi , Kepler-22 b is likely to be a small- 
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Fig. 17. Representation of the liquid water habitable zone around HD 40307 as described by ISelsis et all d2007l) compared to our 
solar system. Note the compactness of the orbital configuration of the five inner planet candidates around HD 40307 compared to 
the planets in the inner solar system. 



scale ve rsion of Neptune rath er than a rocky planet with a solid 
surface dBorucki et all 120121) . However, the precise natures of 
GJ 667C c and HD 40307 g are unknown and they could also 
be scaled-down versions of Neptune-like planets with thick at- 
mospheres and solid cores. Compared to the candidates orbiting 
nearby low-mass stars (GJ 581 d, GJ 667C c, HD 85512 b), HD 
40307 g is not likely to suffer from tidal-locking - which im- 
proves its chances of hosting an Earth-like climate. Like for the 
other HZ candidates, it is not yet possible to determine its phys- 
ical and geochemical properties and a direct-imaging mission 
seems to be the most promising way of obtaining observational 
information on these properties. 

For stable stars, the HARPS data can be readily used 
to detect Keple rian period i cities in RVs with amplitudes be- 
low 1.0 ms _1 (PepeetajJ 1201 ll) . This accuracy can be ad- 
ditionally improved by obtaining RV measurements using 
optimized data-reduction techn iques (e.g. HARPS-TERRA; 
Anglada-Escude & Butlet], 1201 2h and analysing the wavelength 
dependence of the signals. Also, detections should not rely on 
criteria solely based on the analyses of the model residuals - 
this procedure is prone to biases and the weakest signals in the 
data can easily escape detection (compare e.g. th e results of 
lLovis et al 1 1201 laHFeroz et all 1201 lblTuomill2012l) . 

The ability to detect such low-amplitude signals also intro- 
duces a practical problem. The orbits of the planets can no longer 
be readily compared with large numbers of data points using 
phase-folded curves (Fig.[T2b because - while the signals may be 
significant - their existence cannot be so readily verified by eye- 
sight from such figures. Therefore, when discovering low-mass 
planets with the radial velocity method, it is necessary to demon- 
strate the significance of the signals using some other means. 
The detection criterion we describe, namely, that the radial ve- 
locity amplitude is statistically distinguishable from zero, pro- 
vides an easy way of demonstrating this. The distributions of all 
amplitude parameters in Fig.[13]clearly satisfy this criterion and 
show that all signals we report are indeed significantly present 
in the data. 

Additional measurements are needed to verify the existence 
and better constrain the orbits and minimum masses of the three 



new planet candidates reported here. Moreover, given that all 
signals are of planetary origin, detailed dynamical analyses can 
additionally constrain the orbits and masses by excluding unsta- 
ble configurations, and especially, constraining the upper limits 
of the planetary masses that could then be used to determine the 
inclinations of their orbits. In this sense, and given the relatively 
dense dynamical packing of the system, the planet candidates 
around HD 40307 are likely to have masses close to their mini- 
mum values. 

The outer planet candidate (HD 40307 g) has a minimum 
mass of ~ 7 M e and orbits the star within the liquid water HZ. 
Given its priviledged position comfortably within the HZ, it is 
likely to remain inside it for the full orbital cycle even if its ec- 
centricity is not strictly zero. A more detailed characterization 
of this candidate is very unlikely using ground based studies be- 
cause it is very inlikely to transit the star, and a direct imaging 
mission seems the most promising way of learning more about 
its possible atmosphere and life-hosting capabilities. The plane- 
tary system around HD 40307 has an architecture radically dif- 
ferent from that of the solar system (see Figure fTTb. which indi- 
cates that a wide variety of formation histories might allow the 
emergence of roughly Earth-mass objects in the habitable zones 
of stars. 
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Table A.2. HARPS-TERRA Differential radial velocities of HD 40307 from the full spectra (F) and the red part (490-680 nm) of the spectra 
(R). These velocities are in the solar system barycentric reference frame and are corrected for the perspective acceleration effect. The S-index 
measurements in the Mount Wilson system and corresponding uncertainties (photon noise) are also provided. The rightmost column contains 
the measurements of the FWHM of the cross correlation function as provided by the HARPS-DRS. Since no uncertainties for the FWHM are 
provided by the HARPS-DRS, we use the measured standard deviation over all 135 nightly averaged measurements (6.5 m s _1 ) as a measure of 
the uncertainty of each nightly average used in the analysis. 
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Table A.2. Continued. 
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